Machine learning and deep learning approaches typically require extensive simulation or experimental data, which can be difficult to obtain due to the complexity and expense involved in both simulations and experiments. When data is scarce, data-driven methods often struggle to achieve sufficient accuracy. Furthermore, these methods, when faced with noisy real-world observations (such as sensor or experimental data) or incorrectly labeled datasets, tend to make unphysical predictions because there is no information about the underlying physics-based mechanism to validate the results. To address the issue of limited data and the opaque nature of purely data-driven approaches, Physics-Informed Neural Networks (PINNs) have been introduced. PINNs incorporate differential equations and/or other mathematical expressions of physical laws alongside measurement data to the loss function of the neural network [1]. In civil engineering, for example, PINNs can be utilized in hybrid digital twins to simulate structures based on sensor data from physical objects. In a previous study, we have successfully extended the PINN framework for solving solid mechanics and especially contact mechanics problems in 2D [2]. This contribution now presents a further extension of this framework to 3D and time-dependent problems. The PINNs are implemented using a mixed variable formulation, where both displacements and stresses are predicted as network outputs, allowing to formulate the governing equations only as first-order derivatives of the network inputs. Additionally, for time-dependent problems, an interesting approach is further explored, where displacements, stresses, and velocities are all predicted by the neural network. This requires additional coupling operators between the three solution fields. Boundary conditions of the investigated problems are enforced in two ways: as so-called soft constraints and hard constraints. Enforcing soft constraints is achieved by adding specific contributions to the overall loss function while enforcing the hard constraints is accomplished by transforming the network outputs to fulfill the boundary conditions directly. The mixed variable formulation also enables enforcing Neumann boundary conditions and initial conditions as hard constraints. We study our proposed formulation on various examples including the Lamé problem, single-patch test, and the Hertzian contact problem [3].